SYS 2202 final project

Introduction

Related work

Spatiotemporal Characteristics of the COVID-19 Epidemic in the United States

This paper discusses a somewhat similar methodology to the one used in this analysis in that it attempts to use national dashboard data at the county-level to investigate the spatio-temporal dynamics of COVID-19 spread in the United States. However, as the study was conducted earlier in the pandemic (Jan 22 2020 - May 20 2020), it only includes data up to that point. Trends were investigated by geographic region of the US by week, as well.

The paper’s methodology is different. It uses a specially-developed space-time scan statistic constructed to investigate the outbreak of diseases to identify clusters. https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.0020059. This is more in the line of statistical significance testing rather than clustering.

The paper also applied a classification approach (logistic regression) to determine whether a particular county was part of a hotspot or not. The regression coefficients were tested for statistical significance at .05.

The paper also investigated the effects of other regional variables such as population density, distance to nearest airport, % older population, % minority population, and % poverty.

The paper’s analysis suggested that:

• Spatiotemporal trends:
  • 2 geographic clusters (of new cases) were observed in northern and southern Califormia in the 1st 6 epi-weeks.
  • Geographic clusters were observed in Washington, New York, Massachusetts, and Northern Florida. by the 7th epiweek.
  • Following that, Clusters rapidly increased all across the country, with many in the Midwest and South regions.

TODO: INSERT FIGURE OF CLUSTERS HERE

Taking the pulse of COVID-19: a spatiotemporal perspective

This paper is a literature review on the different perspectives that have been taken when performing analyses of COVID-19 time series data. It was received by the publishers in June 24 2020.

• The paper begins by discussing the large number of fundamental uncertainties about the nature of COVID-19 and its spread (e.g reproduction rate, the effects of government policy, the effects of warmer weather)

• The authors propose a so-called “spatiotemporal matrix” for understanding and reconciling COVID-19 data within a larger picture. In the information age, a large volume of “big data” is produced and so a robust conceptual framework is needed to reconcile the information.

  • 5 major sectors of Policy, Environment, Economics, Human Dynamics, Social Implications, and Medical Services.

• Authors discuss the overview of the worldwide spread of COVID-19

  • In particular, they note that in the US:
    • Initially was concentrated on east and west coasts, but spread nearly everywhere through transportation hub regions
    • Feb, Early march saw epicenters in Seattle, LA, and Boston
    • New York became the world epicenter by mid-march
    • By March 17, COVID-19 was present in all states
    • Hot spots emerged in states like California and southern states like Florida, Texas (premature opening) taking the lead

• Additional important aspects are discussed

  • The variation of public policy over time and the developed COVID-19 “Stringency index”
  • Human mobility patterns during COVID-19
  • The initial reduction of emissions during lockdown periods
  • Economic impacts
  • the impact on minority and marginalized communities
  • Possible forecasting techniques.

Progression of COVID‐19 From Urban to Rural Areas in the United States: A Spatiotemporal Analysis of Prevalence Rates

This paper investigated the spread of COVID-19 into rural regions, which is of concern as rural regions typically have less-developed healthcare infrastructure than urban regions. Thus, the authors investigated changes in estimated prevalance rates over time at the county level, comparing between rural-urban status. The data used ranged from April 3, 2020 to April 22, 2020.

Overall, 4 hypotheses were explored: - Investigate county-wide change in prevalence rates by rural-urban status - Identify hotspots based on prevalence rates - Find association of demographic, smoking, disease variables with prevalence rates at the county level - identify counties that show rapid changes in prevalence rates.

• To explore the 1st hypothesis, summary statistics for each county were computed and significance tests were performed to determine if there were significant differences between counties.

• A stochastic model of the spatio-temporal changes in prevalence rates was fitted using MCMC methods and used to denoise the raw data. The estimated prevalence rates were then used to assess the associations between demographic variables.

Time to address the spatiotemporal uncertainties in COVID-19 research: Concerns and challenges

• This paper is a short article detailing some caution that someone performing spatio-temporal analyses of COVID-19 might need to watch out for, especially those performing ecological studies.

• Essentially, temporal and spatial aggregation schemes (i.e, smallest unit of time and smallest unit of area) may bias the data in systematic ways, resulting in analyses that may have different outcomes.

  • Other research has shown that aggregating at large geographic scales may not be granular enough to uncover meaningful associations.

• Ecological studies also suffer from person-level differences which may be a strong confounding variable.

  • Not accounting for the different environments and times that people expose themselves in different environments while trying to account for associations between environmental factors and spread can result in spurious analyses.
  • So called “temporal misalignment”

• The authors of this article advocate for stronger experimental designs, especially individual-level analyses facilitated by anonymized cell-location data.

  • Rather than performing ecological studies, perform more tightly-controlled observational ones.

Main findings:

Analysis design

Question: “How have COVID hotspots changed over time?”

We are interested in understanding the spatial and temporal characteristics of the spread of COVID-19 in the contiguous US, e.g “Where and when do COVID cases and deaths increase or decrease?” To that end, we utilize two approaches to investigate this hypothesis.

• Use choropleth maps (i.e “heatmaps”) to visualize the increase in cases per month in each state. Choropleth maps are intuitive and easy to understand, and can be an effective form of EDA for this problem.
  • We will first create maps at the national level to understand national trends in cases in deaths over time.
  • They will be created with spatial units of states, and temporal units of months in order to understand broad, sweeping trends first before conducting more specific analyses.
  • Then, we will create maps for each region of the contiguous US (e.g Northeast, South, Midwest, etc.) in order to understand trends by region. This analysis is necessary since the national map will not have a scale sensitive enough to capture regional trends.
• Then, use k-means clustering to uncover trends and patterns, comparing them to the choropleth maps to see if they make sense. Along with understanding distinct patterns of spread, k-means may also help us understand if current schemes of segmentation (e.g cases by region) are justified.
  • It may also motivate analyses of some other factor on the spread of COVID-19 at the state level, e.g associations between demographic data, state policy, and the spread of COVID-19.
  • In k-means, we will cluster on time (date, at the granularity of day), space (lat/long, at the granularity of state), and COVID metrics (increase in cases, increase in deaths.)
  • The standard Euclidean distance will be used as the similarity metric, and elbow plots will be used to determine a suitable number of clusters.

Analysis

Part 1: Spatio-temporal heatmapping of cases and deaths across the contiguous US

• Goal: Use animated choropleth maps to investigate how new cases and new deaths change over time intervals.
  • We are interested in broad trends in states. Using a state-level choropleth map is a readily digestible way of seeing how the case and death data change in the US over time.

Data cleaning

library(tidyverse)
library(lubridate)

• The COVID-19 Data Repository by CSSE at JHU was used.
  • Specifically, the time-series csvs were used.
• Retrieve the most recent from the GitHub.

JHUdata_cases <- read_csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv")

JHUdata_deaths <- read_csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv")

• Both datasets have the same structure that is somewhat problematic.
  • The dates are column headings rather than grouped into a single variable.
  • The dataset contains county-level data that will be cleaned when we do state-level plotting.

JHUdata_deaths %>% head()

JHUdata_cases %>% head()

• Perform the same set of cleaning steps:
  • Use pivot_longer with a regex to collapse all the date columns into a date variable.
  • Rename columns to simpler names.
  • Parse the date field from chr to date.
  • Collapse county-level data into state-level data by using group_by and summarize, since we are interested first in state-level trends.
    • This variable is added as cumul for cumulative state-level data.

cases <- JHUdata_cases %>% 
  pivot_longer(matches("[0-9]+/[0-9]+/[0-9]+"),names_to="date",values_to="cumul") %>% 
  rename(c("state"="Province_State")) %>% 
  mutate(date=mdy(date)) %>% 
  group_by(date,state) %>% summarize(cumul=sum(cumul))

deaths <- JHUdata_deaths %>% 
  pivot_longer(matches("[0-9]+/[0-9]+/[0-9]+"),names_to="date",values_to="cumul") %>% 
  rename(c("state"="Province_State")) %>% 
  mutate(date=mdy(date)) %>% 
  group_by(date,state) %>% summarize(cumul=sum(cumul))

• However, we would like to plot new cases over time and the dataset doesn’t contain a “new cases” column. So we need to get a first-order difference of the cumulative time-series.
  • dplyr lag didn’t function properly and returned NAs. So, an in-house function to compute first-order differences was written.
  • first_diff simply iterates through a time series and appends the difference to the correct slot in a vector.
  • add_firstdiff is meant to take a national dataframe and get the 1st order difference for each state.
    • We do this by iterating over every state, filtering the dataframe to get that state, and creating a list of first order differences for each state.
    • Then we combine them all into one dataframe using reduce and full_join.

# Cumulative plot!
# Not what we want.
cases %>% filter(state=="Virginia") %>% 
  ggplot(aes(x=date,y=cumul)) +
  geom_line()

first_diff <- function(arr){
  to_return <- NULL
  for (i in 2:length(arr)){
   to_return[i] <- arr[i] - arr[i-1]
  }
  return(to_return)
}

add_firstdiff <- function(df){
  unique_states <- unique(df$state)
  dfs <- list()
  for (this_state in unique_states){
    thisdf <- df %>% filter(state==this_state)
    thisdf$diff <- first_diff(thisdf$cumul)
    dfs[[length(dfs) + 1]] <- thisdf
  }
  return(dfs %>% reduce(full_join))
}

• The firstdiff functions are applied to the case and deaths dataframes to get the new cases and new deaths per day.
  • They are then merged together with left_join to get a master dataframe.
  • Only the relevant columns (date, state, new cases/new deaths) are selected.
  • The state variable is converted to lower-case to facilitate joining with GeoJSON data later.
  • We filter out all days before Feb 02, 2020 as there were not many cases in January (relative to the largest spikes in the pandemic), so they will not show up very clearly on the heatmap.
    • Also, it makes the time intervals convenient to depict because then we can consider whole months beginning from the 1st day of every month.
  • We make sure to filter out regions that are not US states.
    • The dataset also includes information about US territories (and even a cruise ship!)

cases <- add_firstdiff(cases) %>% 
  select(date,state,diff) %>% 
  rename(c("positiveIncrease"="diff"))
deaths <- add_firstdiff(deaths) %>% 
  select(date,state,diff) %>% 
  rename(c("deathIncrease"="diff"))

data <- left_join(cases,deaths,by=c("date","state")) %>% 
  mutate(state=tolower(state)) %>% 
  filter(date >= date("2020-02-01"))

data <- data %>% filter(state %in% tolower(state.name))

# Our difference function worked!
cases %>% filter(state=="Virginia") %>% 
  ggplot(aes(x=date,y=positiveIncrease)) +
  geom_line()

data %>% ggplot(aes(x=date,y=positiveIncrease)) + geom_line(aes(color=state),alpha=0.3)

• We want to get the dates that our plots will range over.
  • We’ll generate one plot for each month of the pandemic. So use seq to interpolate between the earliest date and the latest date by months.
  • The last date is manually added (it might not be a full month to the last date in the data frame) so that we can also incorporate the most recent information as well.
• We want to generate a set of intervals (ordered pairs) of months, e.g from Jan to Feb, Feb to Mar, Mar to Apr..
  • A simple for-loop takes care of this and stores the info in intervals.

dates <- seq(min(data$date),max(data$date),by="months")
dates <- dates %>% append(max(data$date))

intervals <- list() 
for (i in 2:length(dates)){
  intervals[[i-1]] <- c(dates[i-1],dates[i])
}

• R has built-in geographical data that represents each state as a set of polygons by lat and long.
  • We will use this to plot the US map.

states_map <- map_data("state")

• Here are the functions that plot a map for one specific interval.
  • It works by filtering out all observations in the data in that interval, grouping by state, and then getting the cumulative increase in each state over that interval of time.
  • Then, that data is joined with the map data to produce a plot using ggplot.
  • The gradient scale is then set at an absolute number to ensure that the scale is consistent across all plots produced.
    • This will make sure we can compare maps across time.
  • The map is titled accordingly using the paste function so we know which date interval it is.
  • The plot is cleaned by removing some unnecessary elements.

heatmap_for_givenday_deathIncrease <- function(interval){
  cumul_increase <- data %>% 
    filter(between(date,interval[1],interval[2])) %>% 
    group_by(state) %>% 
    summarize(total=sum(deathIncrease))
  #cumul_increase$statename <-  tolower(state.name[match(cumul_increase$state,state.abb)])
  loc_data <- left_join(states_map,cumul_increase,by=c("region"="state"))
  # Used to see the max value in order to set the gradient. 
  #print(max(cumul_increase$total))
  p <- loc_data %>% 
    ggplot() + 
    geom_polygon(aes(x=long,y=lat,group=group,fill=total),color="grey") +
    scale_fill_gradient(low="white",high="red", limits=c(0,22136)) +
    ggtitle(paste("New deaths: ",interval[1]," to ",interval[2])) +
    theme(axis.title=element_blank(),
          axis.ticks=element_blank(),
          panel.grid=element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank())
  return(p)
}

heatmap_for_givenday_positiveIncrease <- function(interval){
  cumul_increase <- data %>% 
    filter(between(date,interval[1],interval[2])) %>% 
    group_by(state) %>% 
    summarize(total=sum(positiveIncrease))
  #cumul_increase$statename <-  tolower(state.name[match(cumul_increase$state,state.abb)])
  loc_data <- left_join(states_map,cumul_increase,by=c("region"="state"))
  #loc_data <- left_join(states_map,cumul_increase,by=c("region"="statename"))
 
  # Used to see the max value in order to set the gradient. 
  #print(max(cumul_increase$total))
  
  
  p <- loc_data %>% 
    ggplot() + 
    geom_polygon(aes(x=long,y=lat,group=group,fill=total),color="grey") +
    #scale_fill_gradient(low="white",high="red",limits=c(0,1259607)) +
    scale_fill_gradient(low="white",high="red",limits=c(0,1259607)) +
    ggtitle(paste("New cases: ",interval[1]," to ",interval[2])) +
        theme(axis.title=element_blank(),
          axis.ticks=element_blank(),
          panel.grid=element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank())
  return(p)
}

• We can use map from purrr to apply our function readily across our list of intervals, making the process of generating each plot very straightforward.
  • A sample plot is displayed below (there are too many to include in this document).
  • They are exported by running each line again manually.

plots <- map(intervals,heatmap_for_givenday_deathIncrease)
plots <- map(intervals,heatmap_for_givenday_positiveIncrease)
plots[[12]]

• These maps can be exported batch-by-batch using the code below.

for (i in 1:length(plots)){
  ggsave(plot=plots[[i]],path="../project/out",filename=paste("cases",i,".png",sep=""),width=5,height=4,device="png")
}

• ffmpeg was then used to join the images together into a video.
  • The command needed to be adjusted slightly from batch to batch.

ffmpeg -r 1 -f image2 -s 1920x1080 -i ../deaths/project/out/cases/cases%d.png -vcodec libx264 -crf 25 -pix_fmt yuv420p test.mp4

• The videos cannot be embedded into a PDF,

• This whole procedure can also be replicated for each region of the US, giving us a better idea of regional trends as well.

  • A csv containing the region names for each state is imported and joined with the data.

us_regions <- read_csv("https://raw.githubusercontent.com/cphalpert/census-regions/master/us%20census%20bureau%20regions%20and%20divisions.csv")

data <- left_join(cases,deaths) %>% 
  filter(date >= date("2020-02-01"))

data <- left_join(data,us_regions,by=c("state"="State"))

• The functions are adjusted slightly to allow filtering by region.
  • An upbound argument is added to make it easier to set the absolute scale of the gradient.

heatmap_for_givenday_positiveIncrease <- function(interval,region,upbound){
  cumul_increase <- data %>% 
    filter(Region == region) %>% 
    filter(between(date,interval[1],interval[2])) %>% 
    group_by(state) %>% #added when JHU data used
    mutate(state=tolower(state)) %>% 
    summarize(total=sum(positiveIncrease))
  #cumul_increase$statename <-  tolower(state.name[match(cumul_increase$state,state.abb)])
  #loc_data <- inner_join(states_map,cumul_increase,by=c("region"="statename"))
  loc_data <- inner_join(states_map,cumul_increase,by=c("region"="state"))
  
  print(max(cumul_increase$total))
  
  p <- loc_data %>% 
    ggplot() + 
    geom_polygon(aes(x=long,y=lat,group=group,fill=total),color="grey") +
    scale_fill_gradient(low="white",high="red",limits=c(0,upbound)) +
    ggtitle(paste(region,"- New cases: ",interval[1]," to ",interval[2])) +
        theme(axis.title=element_blank(),
          axis.ticks=element_blank(),
          panel.grid=element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank())
  return(p)
}

heatmap_for_givenday_deathIncrease <- function(interval,region,upbound){
  cumul_increase <- data %>% 
    filter(Region == region) %>% 
    filter(between(date,interval[1],interval[2])) %>% 
    group_by(state) %>% 
    mutate(state=tolower(state)) %>% 
    summarize(total=sum(deathIncrease))
  #cumul_increase$statename <-  tolower(state.name[match(cumul_increase$state,state.abb)])
  #loc_data <- inner_join(states_map,cumul_increase,by=c("region"="statename"))
  loc_data <- inner_join(states_map,cumul_increase,by=c("region"="state"))
  
  print(max(cumul_increase$total))
  
  p <- loc_data %>% 
    ggplot() + 
    geom_polygon(aes(x=long,y=lat,group=group,fill=total),color="grey") +
    scale_fill_gradient(low="white",high="red",limits=c(0,upbound)) +
    ggtitle(paste(region,"- New deaths: ",interval[1]," to ",interval[2])) +
        theme(axis.title=element_blank(),
          axis.ticks=element_blank(),
          panel.grid=element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank())
  return(p)
}

• The maps are once again generated in a very similar procedure.
  • Below is an example regional map.

plots <- map(intervals,heatmap_for_givenday_positiveIncrease,"South",639843)
plots <- map(intervals,heatmap_for_givenday_positiveIncrease,"West",1259607)
plots <- map(intervals,heatmap_for_givenday_positiveIncrease,"Northeast",450056)
plots <- map(intervals,heatmap_for_givenday_positiveIncrease,"Midwest",322299)

plots <- map(intervals,heatmap_for_givenday_deathIncrease,"South",10159)
plots <- map(intervals,heatmap_for_givenday_deathIncrease,"West",15922)
plots <- map(intervals,heatmap_for_givenday_deathIncrease,"Northeast",22136)
plots <- map(intervals,heatmap_for_givenday_deathIncrease,"Midwest",6538)

plots[[12]]

• A similar code and bash command is used to generate animations.

for (i in 1:length(plots)){
  ggsave(plot=plots[[i]],path="../project/out",filename=paste("northeast_cases",i,".png",sep=""),width=5,height=4,device="png")
}

• National cases:

• National deaths:

• Midwest cases

• Midwest deaths

• Northeast cases

• Northeast deaths

• South cases

• South deaths

• West cases

• West deaths

Part 2: Spatio-temporal clustering of cases and deaths across the contiguous US at the state level

• The animated maps are a great way for us to quickly see broad trends in cases and deaths across the US. We can also try using k-means clustering as a data-mining approach to see if we can uncover additional trends and patterns by clustering.
  • This might help us understand if there are some “phases” or “clusters” of the pandemic in the US, not only just in space but also across time.

# Make analysis reproducible
set.seed(2)
library(tidyverse)
library(lubridate)
library(cluster)
library(rvest)
library(patchwork)

• Perform roughly the same data cleaning steps as before.
  • Be sure to omit Alaska and Hawaii, as those are non-contiguous and will interfere with the cluster analysis.

JHUdata_cases <- read_csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv")

JHUdata_deaths <- read_csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv")

first_diff <- function(arr){
  to_return <- NULL
  for (i in 2:length(arr)){
   to_return[i] <- arr[i] - arr[i-1]
  }
  return(to_return)
}

add_firstdiff <- function(df){
  unique_states <- unique(df$state)
  dfs <- list()
  for (this_state in unique_states){
    thisdf <- df %>% filter(state==this_state)
    thisdf$diff <- first_diff(thisdf$cumul)
    dfs[[length(dfs) + 1]] <- thisdf
  }
  return(dfs %>% reduce(full_join))
}

cases <- JHUdata_cases %>% pivot_longer(matches("[0-9]+/[0-9]+/[0-9]+"),names_to="date",values_to="cumul") %>% 
  rename(c("state"="Province_State")) %>% 
  mutate(date=mdy(date)) %>% 
  group_by(date,state) %>% summarize(cumul=sum(cumul))

deaths <- JHUdata_deaths %>% pivot_longer(matches("[0-9]+/[0-9]+/[0-9]+"),names_to="date",values_to="cumul") %>% 
  rename(c("state"="Province_State")) %>% 
  mutate(date=mdy(date)) %>% 
  group_by(date,state) %>% summarize(cumul=sum(cumul))


cases <- add_firstdiff(cases) %>% select(date,state,diff) %>% rename(c("positiveIncrease"="diff"))
deaths <- add_firstdiff(deaths) %>% select(date,state,diff) %>% rename(c("deathIncrease"="diff"))

data <- left_join(cases,deaths) %>% 
  ungroup() %>%
    filter(state %in% state.name) %>% 
    filter(state != "Alaska") %>%
    filter(state != "Hawaii")

• Check the dataframe and produce the same plot to make sure everything went smoothly.

data %>% head()

data %>% ggplot(aes(x=date,y=positiveIncrease)) + geom_line(aes(color=state),alpha=0.3)

• Although the JHU dataframe does contain location data, the locations provided are at the county level and we are looking to do a clustering at the state level.
  • Full disclosure, I attempted to do a county-level clustering. However, I was unable to get valid county-level 1st differences, as due to retroactive corrections of death and case data, extremely large negative values were observed.
  • Thus, since the cleaning steps for the state level had already been performed, I opted to continue with a state-level analysis.
• So, the state lat-long data must be added to include a spatial component in the clustering. The location data is scraped from the web and joined onto the data.
  • The locations provided appear to be roughly central in each state, most likely the capital cities.

content <- read_html("https://inkplant.com/code/state-latitudes-longitudes")
states_latlong <- html_table(content,header=T)

states_latlong <- as_tibble(data.frame(states_latlong))
states_latlong <- states_latlong %>% rename(c("state"="State","lat"="Latitude","long"="Longitude"))

# Check how the states distribute
# Looks like the capital cities of each state, that's fine
states_map <- map_data("state")
ggplot(states_map, aes(long, lat)) +
  geom_polygon(aes(group = group),fill = "white", colour = "black") +
  geom_point(data=states_latlong)

data <- data %>% left_join(states_latlong)
data %>% head()

• Now, for clustering. We are going to use the default Euclidean distance as our similarity metric, as we are going to cluster on date/time, new cases, new deaths, and latitude and longitude. Euclidean distance makes sense for these continuous values.

• However, we need to scale or normalize each feature before we use k-means. This is done so that the differing units of the features won’t affect the similarity between 2 observations. Thus, min-max normalization was used to scale the features.

Min-max: To scale a vector \(X\), \(x_{i_{scaled}} = \frac{x_i - min(X)}{max(X) - min(X)}\) for each \(x_i\) in \(X\).

# Function to normalize a vector
norm <- function(my_vec){
  my_max <- max(my_vec)
  my_min <- min(my_vec)
  if (my_max-my_min ==0){
    print("err: divide by zero")
  }
  return((my_vec - my_min)/(my_max-my_min))
}

• There is a fairly strong correlation between deaths and cases. Thus, we probably won’t want to include both at the same time in the clustering analysis. We can run individual analyses with both and see what patterns we observe.

temp <- data %>% filter(!is.na(positiveIncrease) & !is.na(deathIncrease)) 
cor(temp$positiveIncrease,temp$deathIncrease)

## [1] 0.6295024

# Prepare the data for clustering, i.e
# remove NAs,
# and normalize features
cluster_deaths <- data %>% 
filter(
  !is.na(positiveIncrease) & 
  !is.na(deathIncrease) &
  !is.na(lat) & 
  !is.na(long)) %>%
mutate(date=as.integer(date)) %>%
mutate(date=norm(date)) %>%
mutate(
  positiveIncrease=norm(positiveIncrease),
  deathIncrease=norm(deathIncrease),
  lat=norm(lat),
  long=norm(long)
) %>% select(-state,-positiveIncrease)

cluster_cases <- data %>% 
filter(
  !is.na(positiveIncrease) & 
  !is.na(deathIncrease) &
  !is.na(lat) & 
  !is.na(long)) %>%
mutate(date=as.integer(date)) %>%
mutate(date=norm(date)) %>%
mutate(
  positiveIncrease=norm(positiveIncrease),
  deathIncrease=norm(deathIncrease),
  lat=norm(lat),
  long=norm(long)
) %>% select(-state,-deathIncrease)


cluster_deaths %>% head()

cluster_cases %>% head()

• A question that always needs to be answered during k-means is the choice of \(k\). Here, we use the elbow plot approach to get an idea of what k to choose.

# Generate an elbow plot to determine good value of k
plot_elbow <- function(data,up_to){
  elbow <- unlist(map(1:up_to, function(k){
    clusters <- kmeans(data,k,50)
    return (clusters$tot.withinss)
  }))
  data.frame(k=1:up_to,elbow) %>% ggplot(aes(x=k,y=elbow)) +
    geom_line() + geom_point(size=5) +
    scale_x_continuous(breaks=1:20)
}

• The elbow plots produced are overall pretty similar and are a bit unclear (no “dramatic” elbows.)
  • We see a slight elbow at \(k=3\) for cases.
  • For deaths, there is another slight elbow at \(k=6\).
  • \(k=3\) in deaths isn’t really too much of an elbow.
• Despite the lack of a definitive plot, we still have an idea of what k to use. It’s probably best to use no more than \(k=6\), so we’ll try both \(k=3\) and \(k=6\) for additional segmentation.
  • Even though the elbow plots are not particularly uniform across the 2, I chose to use identical values of k to facilitate comparison between the two when clustering.

plot_elbow(cluster_deaths,15) + ggtitle("Cases") + 
  plot_elbow(cluster_cases,15) + ggtitle("Deaths")

clusters_deaths_3 <- kmeans(cluster_deaths,centers=3,50)
clusters_cases_3 <- kmeans(cluster_cases,centers=3,50)

clusters_deaths_6 <- kmeans(cluster_deaths,centers=6,50)
clusters_cases_6 <- kmeans(cluster_cases,centers=6,50)

# Add the cluster IDs to the data.
add_clust_nums <- function(clust_out){
  data_w_clust <- data %>% 
  filter(
    !is.na(positiveIncrease) & 
    !is.na(deathIncrease) &
    !is.na(lat) & 
    !is.na(long))
  data_w_clust$clust_num <- clust_out$cluster
  return(data_w_clust)
}

wclust_deaths_3 <- add_clust_nums(clusters_deaths_3)
wclust_deaths_6 <- add_clust_nums(clusters_deaths_6)
wclust_cases_3<- add_clust_nums(clusters_cases_3)
wclust_cases_6 <- add_clust_nums(clusters_cases_6)

• Visualizing the clusters
  • We can create time-series plots of each cluster by aggregating deaths and cases in that cluster, in order to see the differences between clusters over time.
  • We also want to explore within-cluster variation in new cases and deaths over time, to determine if the clustering algorithm’s choices are reasonable and make sense.
  • We also want to create maps of the clusters to see if there is a shift in location over time.

days <- seq(min(data$date),max(data$date),by="months")
days[length(days)+1] <- max(data$date)

colors <- c("1"="#E69F00", "2" = "#56B4E9", "3" = "#009E73", "4" = "#F0E442", "5" ="#0072B2", "6" ="#D55E00")


plot_clust_timeseries <- function(w_clust_i,name){
  loc_data <-  w_clust_i %>% group_by(date,clust_num) %>%
    summarize(cases=sum(positiveIncrease), deaths=sum(deathIncrease)) 
  
  p1 <- loc_data %>% ggplot(aes(x=date,y=cases)) + 
    geom_smooth(aes(color=factor(clust_num))) + ggtitle(paste("Time series of cases by cluster:",name)) +
    scale_x_continuous(breaks=days[c(T,F)]) +
    scale_color_manual(values=colors)
    
  p2 <- loc_data %>% ggplot(aes(x=date,y=deaths)) + 
    geom_smooth(aes(color=factor(clust_num))) + ggtitle(paste("Time series of deaths by cluster:",name)) +
    scale_x_continuous(breaks=days[c(T,F)]) +
    scale_color_manual(values=colors)
  return(p1 / p2)
}

• Cluster trends

plot_clust_timeseries(wclust_deaths_3,"deaths_3")

plot_clust_timeseries(wclust_deaths_6,"deaths_6")

plot_clust_timeseries(wclust_cases_3,"cases_3")

plot_clust_timeseries(wclust_cases_6,"cases_6")

- Within-cluster trends

wclust_cases_3 %>%  ggplot(aes(x=date,y=positiveIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in positive cases: cases_3")

wclust_cases_3 %>% ggplot(aes(x=date,y=deathIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in deaths: cases_3")

wclust_cases_6 %>% ggplot(aes(x=date,y=positiveIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in positive cases: cases_6")

wclust_cases_6 %>% ggplot(aes(x=date,y=deathIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in deaths: cases_6")

wclust_deaths_3 %>% ggplot(aes(x=date,y=positiveIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in positive cases: deaths_3")

wclust_deaths_3 %>% ggplot(aes(x=date,y=deathIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in deaths: deaths_3")

wclust_deaths_6 %>% ggplot(aes(x=date,y=positiveIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in positive cases: deaths_6")

wclust_deaths_6 %>% ggplot(aes(x=date,y=deathIncrease,color=state,alpha=0.2)) + geom_smooth() + facet_wrap(~clust_num) + ggtitle("Within-cluster time-series, increase in deaths: deaths_6")

states_map <- map_data("state")
cluster_plot_for_day <- function(day,data_w_clust,name){
  temp <- data_w_clust %>% filter(date==day) %>%
    mutate(state=tolower(state)) %>% select(-lat,-long)
  loc_data <- left_join(states_map,temp,by=c("region"="state"))

  loc_data %>% 
    ggplot() + 
    geom_polygon(aes(x=long,y=lat,group=group,fill=factor(clust_num)),color="grey") +
    scale_fill_manual(values=colors) +
    ggtitle(paste("Location of clusters",day,name,sep=" - ")) +
    theme(axis.title=element_blank(),
          axis.ticks=element_blank(),
          panel.grid=element_blank(),
          axis.text.x = element_blank(),
          axis.text.y = element_blank())
}

• In order to visualize the clusters, a couple of representative dates are selected, and we can create a choropleth map by coloring each state according to the cluster it belongs to.

cluster_plot_for_day("2020-07-01",wclust_cases_3,"cases_3")

cluster_plot_for_day("2021-01-01",wclust_cases_3,"cases_3")

cluster_plot_for_day("2020-07-01",wclust_cases_6,"cases_6")

cluster_plot_for_day("2021-01-01",wclust_cases_6,"cases_6")

cluster_plot_for_day("2020-07-01",wclust_deaths_3,"cases_3")

cluster_plot_for_day("2021-01-01",wclust_deaths_3,"deaths_3")

cluster_plot_for_day("2020-07-01",wclust_deaths_6,"deaths_6")

cluster_plot_for_day("2021-01-01",wclust_deaths_6,"deaths_6")

Analysis and discussion

• Several key trends were observed in the animated choropleth maps.

At the national level, the spread of COVID is indeed generally characterized a spread from East coast/Northeast states and West coast states through the entire country. Large spikes of cases continue to be concentrated in larger, populous states like New York and California, although strong increases are also observed in Southern and Midwestern states.

1. The state of New York was a major COVID hotspot during the earlier periods of the pandemic, exhibiting both a high case and death load. New cases and new deaths reached their highest point in this wave from 2020-04-01 to 2020-05-01. The observed death load in this month was one of the highest ever new deaths during the entire pandemic in the contiguous US thus far. At a glance, it seems to be disproportionately high (in terms of the confirmed death/case ratio) compared to the rest of the pandemic.

• Other states adjacent to New York (New Jersey, Pennsylvania, Massachusetts) also see rises in deaths.
• Midwestern states such as Illinois, Michigan, Ohio also see rises in deaths.
• California and Texas

2. Both cases and deaths began to rise in California, Texas, and Florida around 2020-08 to 2020-09. By this point, almost all states are seeing new cases. There is a small decrease in the new number of cases in these 3 states from 2020-09 to 2020-11.

3. There is a dramatic, almost uniform increase in cases across the entire nation from 2020-11 through 2021-03, Some of the largest increases (Roughly in descending order), are in California, Texas, New York, Florida, Illinois, Ohio, Pennsylvania. Also of note are Arizona, Georgia, North Carolina, and Tennessee. Almost all of these increases are accompanied with roughly proportional rises in deaths.

• California reaches one of the largest monthly increases in cases and deaths from 2021-01 to 2021-02, followed by similarly large increases in (roughly descending order) in Texas, Florida, and New York.

4. New cases begin to decrease from 03-2021 to 04-2021. At this point, the data ends, so it is not possible to say what much about what happens in this time frame.

Region by region,

Midwest: - Illinois, Michigan see the worst spike in deaths from 2020-05 to 2020-06. - Coinciding with the nationwide increase from 2020-11 to 2021-01, Illinois, Ohio, Michigan, and Indiana see large rises in cases and deaths. The rest of the Midwest (Excluding North Dakota, South Dakota, Nebraska) all see a similar rise in cases to each other (that is less than that of Ohio and Illinois). - Cases decrease in almost all Midwestern states from 03-2021 to 04-2021, except for Michigan.

Northeast: - Corresponding to the national map, New York sees considerable rises in cases and deaths. - when New York rises in cases, it is typically accompanied by a rise in its adjacent states, like Pennsylvania, New Jersey, and Massachusetts. - Compared to the rest of the Northeast, Maine, New Hampshire, and Vermont do not see as large (raw count) of a rise in cases.

South: - Corresponding with the trends in the national map, Texas and Florida dominate much of the increases in the South. - The spread throughout the rest of the South is relatively uniform in the sense that during periods of increased spread, most Southern states see similar magnitudes of increases. - However, during the 2020-12 to 2021-03 wave, cases and deaths seemed to be more concentrated in Georgia, Tennessee, North Carolina, and Virginia.

West: - California more or less dominates most of the case spread in the West over the entire range of available time. - Arizona does also at times see considerable rises in deaths and cases as well. It is probably the “2nd hotspot” in the West.

Cluster analysis:

• Overall, clustering by deaths or by cases seemed to have resulted in almost identical clusters, thus justifying the decision to cluster on one or the other.

• The spatial distribution of the clusters supports the practice of dividing the US into its political geographical regions. The 6-cluster runs have clusters that correspond roughly to the Northeast, the South, and the West.

  • The 3 cluster runs seem to divide right down the middle of the US, likely due to the differences between the East and West coasts in terms of spread dynamics.

Investigating the “time series by cluster” plots: - All clustering runs seemed to have created the same temporal slice at around 2020-08. We see a peak before this date and a peak after, indicating that this clustering algorithm probably thinks this is about the right date to delineate the “first wave” and the “second wave.” - This lines up with what we saw in the heatmaps, as we saw the onset of a large increase in cases around that time.

• The plots also all show an interesting characteristic:
  • Overall, the cases before 2020-08 compared to the cases after 2020-08 are lower.
  • But, the deaths before 2020-08 compared to the deaths after are comparable.
  • This may line up with what we saw in the maps: New York having a high amount of deaths despite having a lower amount of cases relative to later in the pandemic.
• The 6 cluster runs add a little more granularity to the trend in the 1st wave:

  • When checking cases_6, the cases of cluster 1 go down while the cases of cluster 5 go up. These correspond, respectively, to the Northeast and the South.

  • This checks out with what we observed in the maps - we saw a strong initial wave in the Northeast, followed by increases in the South shortly after.

  • After this, though, the clusters for the 2nd wave all seem to uniformly go up, indicating that the 2nd wave was a nationwide event.

Investigating the within-cluster plots, - the 3-cluster runs seem to have middling results with respect to the time series: in cluster 1, for example, California appears to have had such a large increase in cases that it essentially gets its own cluster. However, the trends of all states within each cluster seem to generally correspond to one another.

• Again, in the 6 cluster runs, California is very much a strong outlier within its cluster.

• We see an interesting trend insofar as the “observed similarity” of within cluster trends - for those clusters that are temporally situated earlier in the pandemic, the within cluster variance seems to be larger (the series values have larger vertical distances).

• With the exception of the cluster containing California, the other clusters later in the pandemic seem to have less within cluster variance.

• This may support the hypothesis that over time, pandemic dynamics within the contiguous US become more homogeneous across the entire nation.

The trends observed in the cluster analysis correspond well to those observed using heatmapping, providing extra support for the observed patterns.

These patterns are an excellent stepping stone to subsequent analyses of state-level COVID data. Since this analysis is fundamentally exploratory, it has generated several questions that would be of interest to test in a confirmatory manner. - e.g Since this analysis has suggested that grouping states by their geographical region is indeed logical, we may wish to test that hypothesis with a statistical test (i.e investigate if there is indeed a group difference between regions) -

Ranking regions by record high monthly case volume in ascending order,

region	monthly record high case increase
Midwest	312,326
Northeast	404,981
South	639,843
West	1,249,607

Ranking top 10 states by record high monthly case volume,

Investigating hotspots by region,

• South
  • Florida and Texas are hotspots of the South
• West
  • California and Arizona
• Northeast
  • New York and Pennsylvania
• Midwest
  • Ohio, Illinois

National trends,

Regional trends,

The cases map:

• Most states have no data in the beginning of the pandemic. The 1st cases appear in Washington.
• The region that has the most heat in the 2020-04-13 to 2020-05-13 period is New York.
  • Most other states are not seeing proportionally large rises in cases yet.
  • California is tinged red, as well as some other states in the Northeast.
• Over the next sets of months, we see increases in cases.
  • California becomes more red.
  • Texas becomes more red.
  • Illinois becomes more red.
• In the 2020-06-13 to 2020-07-13 range,
  • California, Texas, and Florida are more red.
    • Cases increase in these states especially during the next month as well.
  • The Southeast is one of the darker regions in the map
  • New York is less red again.
• 2020-08-13 to 2020-09-13
  • The increase in new cases decreases overall.
  • California, Texas, and Florida are still the reddest states.
  • The entire east coast appears to be a light red, signifying that COVID spread is probably well concentrated in areas with higher population.
• 2020-10-13 to 2020-11-13
  • Cases begin to increase across the entire US. This trend worsens through the next month.
    • California, Texas, Florida, and Ohio are some of the more red states.
  • Texas and Illinois are the reddest states.
• 2020-12-13 to 2020-01-13
  • California is especially red. This is the highest new case count in the entire series of plots.
  • Texas, Florida, and new York are some of the other hot spots.
  • The trend eases off in the next months, however, this is where the data set ends so it may be due to a lack of information.

Deaths map: - 2020-03-13 to 2020-04-13 - New York has the worst count of deaths by far. Other states in the Northeast, California are seeing rises as well. - Deaths continue to worsen in New York through the next month.

• 2020-05-13 to 2020-06-13
  • Deaths are reduced from the last month, however the Northeast, California, Illinois are still standing out.
  • Overall deaths have increased in other states.
• 2020-07-13 to 2020-08-13
  • Deaths begin to increase in the South, especially Florida and Texas.
  • Arizona, California, Nevada see increases in deaths.

Overall at the national level, - New York was a hotspot of COVID activity during the earlier period of the pandemic in the US (from March 2020 to June 2020.) - Relatively speaking, deaths were very high during this period. - New cases and deaths slow from 05 2020 to 06-2020. - Then, in 07-2020 to 08-2020, New cases and deaths begin to increase especially in - California, - Texas, - and Florida. - This continues for a bit until 10-2020, where new cases and deaths begin to slow again.

• From 10-2020 to 1-2021, we see especially large nationwide increases.
  • California sees especially large increases in cases and deaths, followed by Texas and Florida.
  • Ohio, illinois, New York, and most of the Northeast and South also see increases.
• Cases and deaths taper on the last couple available months of the dataset (but may be due to lack of data availability.)

Qualitatively speaking, the deaths and the cases charts coincide quite well. - It is known that there is a “lag time” between increases in cases and increases in deaths. This may be due to the time required to report the data as well as the time required for cases to develop. - This is especially evident from December 2020 to January 2021. The cases hit their peak in California in December, but the deaths do not peak until the month after, even though new cases in California go down during that period.